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Abstract. In the present article, we derive an explicit expression for 
the truncated mean and variance for the multivariate normal distribu- 
tion with arbitrary rectangular double truncation. We use the moment 
generating approach of Tallis (1961) and extend it to general uw, 4 and 
all combinations of truncation. As part of the solution, we also give 
a formula for the bivariate marginal density of truncated multinormal 
variates. We also prove an invariance property of some elements of the 
inverse covariance after truncation. Computer algorithms for computing 
the truncated mean, variance and the bivariate marginal probabilities 
for doubly truncated multivariate normal variates have been written in 
R and are presented along with three examples. 
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1 Introduction 


The multivariate normal distribution arises frequently and has a wide range 
of applications in fields such as multivariate regression, Bayesian statistics and 
the analysis of Brownian motion. One motivation to deal with moments of the 
truncated multivariate normal distribution comes from the analysis of special fi- 
nancial derivatives (“auto-callables” or “Expresszertifikate”) in Germany. These 
products can expire early depending on some restrictions of the underlying tra- 
jectory, if the underlying value is above or below certain call levels. In the frame- 
work of Brownian motion the finite-dimensional distributions for log returns at 
any d points in time are multivariate normal. When some of the multinormal 
variates X = (#1,...,2a)' ~ N(p,™) are subject to inequality constraints (e.g. 
a; < x; < b;), this results in truncated multivariate normal distributions. 
Several types of truncations and their moment calculation have been de- 
scribed so far, for example the one-sided rectangular truncation x > a (Tallis, 
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1961), the rather unusual elliptical and radial truncations a < x2’Ra < b (Tallis, 
1963) and the plane truncation Ca > p (Tallis, 1965). Linear constraints like 
a < Ca < bcan often be reduced to rectangular truncation by transformation 
of the variables (in case of a full rank matrix C: a* = C-'a<a<C lb=3*), 
which makes the double rectangular truncation a < a < b especially important. 


The existing works on moment calculations differ in the number of variables 
they consider (univariate, bivariate, multivariate) and the types of rectangular 
truncation they allow (single vs. double truncation). Single or one-sided trunca- 
tion can be either from above (a < a) or below (a > a), but only on one side 
for all variables, whereas double truncation a < x < b can have both lower and 
upper truncation points. Other distinguishing features of previous works are fur- 
ther limitations or restrictions they impose on the type of distribution (e.g. zero 
mean) and the methods they use to derive the results (e.g. direct integration or 
moment-generating function). 


Rosenbaum (1961) gave an explicit formula for the moments of the bivariate 
case with single truncation from below in both variables by direct integration. 
His results for the bivariate normal distribution have been extended by Shah 
and Parikh (1964), Regier and Hamdan (1971) and Muthén (1990) to double 
truncation. 


For the multivariate case, Tallis (1961) derived an explicit expression for the 
first two moments in case of a singly truncated multivariate normal density with 
zero mean vector and the correlation matrix R using the moment generating 
function. Amemiya (1974) and Lee (1979) extended the Tallis (1961) derivation 
to a general covariance matrix & and also evaluated the relationship between 
the first two moments. Gupta and Tracy (1976) and Lee (1983) gave very simple 
recursive relationships between moments of any order for the doubly truncated 
case. However, except for the mean, there are fewer equations than parameters. 
Therefore, these recurrent conditions do not uniquely identify moments of order 
> 2 and are not sufficient for the computation of the variance and other higher 
order moments. 


Table 1 summarizes our survey of existing publications dealing with the com- 
putation of truncated moments and their limitations. Even though the rectangu- 
lar truncation a < a < bcan be found in many situations, no explicit moment 
formulas for the truncated mean and variance in the general multivariate case 
of double truncation from below and/or above have been presented so far in the 
literature and are readily apparent. The contribution of this paper is to derive 
these formulas for the first two truncated moments and to extend and generalize 
existing results on moment calculations from especially Tallis (1961), Lee (1983), 
Leppard and Tallis (1989), and Muthén (1990). Besides, we also refer Kan and 
Robotti (2017) and Arismendi (2013) for the moment computation of folded and 
truncated multivariate normal distribution. However, for moments computation 
for skewed and extended skew normal distribution, we refer Kan and Robotti 
(2017) and Arellano-Valle and Genton (2005). In the sequel, we also make a note 
on the existing R package ” MomTrunc” (see Galarza C.E. & V.H., 2021) for nu- 
merical computation of moments of folded and truncated multivariate normal 
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Table 1. Survey of previous works on the moments for the truncated multivariate 
normal distribution 








Author #Variates Truncation Focus 

Rosenbaum bivariate single moments for bivariate normal vari- 

(1961) ates with single truncation, b; < 
yi <0, b2 < y2< co 

Tallis (1961) multivariate single moments for multivariate normal 
variates with single truncation from 
below 

Shah and Parikh _ bivariate double recurrence relations between mo- 

(1964) ments 

Regier and bivariate double an explicit formula only for the case 

Hamdan (1971) of truncation from below at the 
same point in both variables 

Amemiya (1974) multivariate single relationship between first and sec- 
ond moments 

Gupta and Tracy multivariate double recurrence relations between mo- 

(1976) ments 

Lee (1979) multivariate single recurrence relations between mo- 
ments 

Lee (1983) multivariate double recurrence relations between mo- 
ments 

Leppard and multivariate single moments for multivariate normal 

Tallis (1989) distribution with single truncation 

Muthén (1990) bivariate double moments for bivariate normal dis- 


tribution with double truncation, 
br < yt < a1, b2 < yo < ag 
Manjunath and multivariate double moments for multivariate normal 
Wilhelm distribution with double truncation 
in all variablesa<a<b 





and Student’s t-distribution. However, the aforementioned package also suggests 
the ”tmvtnorm” package (e.g., Wilhelm & Manjunath, 2012), which is solely 
based on the results presented in this note. Finally, we also refer Genz (1992) 
for the numerical computation for the multivariate normal probabilities. 


The rest of this paper is organized as follows. Section 2 presents the moment 
generating function (m.g.f) for the doubly truncated multivariate normal case. 
In Section 3, we derive the first and second moments by differentiating the m.g.f. 
These results are completed in Section 4 by giving a formula for computing the 
bivariate marginal density. In Section 5, we present two numerical examples and 
compare our results with simulation results. Section 6 links our results to the 
theory of graphical models and derives some properties of the inverse covariance 
matrix. Finally, Section 7 summarizes our results and gives an outlook for further 
research. 


20 B. G. Manjunath & S. Wilhelm 


2 Moment Generating Function 


The d-dimensional normal density with location parameter vector yp € R@ and 
non-singular covariance matrix & is given by 





Pu.=(x) = sprTE ex ; (w — ps)’ Zo" (a - »)} , @eER*. (1) 


The pertaining distribution function is denoted by ®,,5(x). Correspondingly, 
the multivariate truncated normal density, truncated at a and b, in R®, is defined 
as 


Pp.E(L) f 
—§~—,,, fora<a<b, 
Popu(#) = 4 P{asX<b} (2) 
0, otherwise. 


Denote a = P {a < X < b} as the fraction after truncation. 

The moment generating function of a d-dimensional truncated random vari- 
able X, truncated at a and b, in R%, having the density f(a) is defined as the 
d-fold integral of the form 


b 
m(t) = E (c! *) = | et © #(x) dex. 
a 
Therefore, the m.g.f for the density in (2) is 


1 


m(t) = a(2m)4/2| 51/2 


ie exp {-3 [(a@ — ps)’ D7! (@ — p) — 2t'ar] \ dx. (3) 


a 


In the following, the moments are first derived for the special case = O. Later, 
the results will be generalized to all yw by applying a location transformation. 
Now, consider only the exponent term in 3 for the case uw = 0. Then we have 


— [2’=" ta — 2t'a] 
which can also be written as 
st Bt 5 [(@-9 D4 (@-8)], 
where € = Ht. 


Consequently, the m.g.f of the rectangularly doubly truncated multivariate 
normal is 


ef b 
m() = tes | ew {se Eee hde, 


where T = stXt. 
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The above equation can be further reduced to 


e! b-§ 1 a 
m(t) = ACS ee ae exp {5a x x} dx. (5) 
For notational convenience, we write equation 5 as 
m(t) = eT Gay (6) 


where 


@ ! [ve 1 le dx 
a = oe Xp 4—= : 
2 a@mPEP? Jae P12 


3 First and Second Moment Calculation 


In this section, we derive the first and second moments of the rectangularly 
doubly truncated multivariate normal density by differentiating the m.g.f.. 

Consequently, by taking the partial derivative of (6) with respect to t; we 
have 








(7) 


In the above equation the only essential terms that will be simplified are 


OeT a : 
at; =e a Oiktk 
k=1 





and 





ey rs) in ie 
= a Pay (x)dxg...dxi, (8) 
Ot; Ot; at a* 

where a = a; — 4 Oint, and bf = b; — Sra J;,~t,. Subsequently, (8) is 

d 


ox = So oin (Fe(aj,) — Fe (8), (9) 
s k=1 





where 


Fj(«) = 
by bit poi gs ba 

/ sai : a Pas (1, --)Vi-1, ©, Vi4i, -Ca)dvg...dxj41duj_1...d21. 
ay aii VOiyy a4 


(10) 
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Note that at t, = 0, for all & = 1,2,...,d, we have a¥ = a; and bF = bj. 
Therefore, F;(a) will be the i-th marginal density. An especially convenient way 
of computing these one-dimensional marginals is given in Cartinhour (1990). 

From (7) — (9) for & = 1,2,...,d all t, = 0. Hence, the first moment is 





B(x) = OY i=) Benue =m). (11) 


k=1 


Now, by taking the partial derivative of (7) with respect to t;, we have 


Pm(t)  pP@an , and Pel. de? Aas, 
tat; | Otjdt; | Ot Ot; | °7OtjO% | Ot Ot; 





(12) 


The essential terms for simplification are 


O7eT 
—— SO 
Ot; Ot; a 


and clearly, the partial derivative of 9 with respect to t; gives 


PPOs. OF; (ax) e OF, (bi) 
aia (a) De (eee) a 


k=1 





In the above equation, merely consider the partial derivative of the marginal 
density Fi,(a;) with respect to t;. With further simplification, it reduces to 


OF, (a bea beat 
ay re X15 0+ UK 1, 0%, 2k41,- Lq)da_ k 
me = 55 


Ohya 


= tins 


Ok,k 
Tk,q?j,k * * * L* 
+35 (<i = ae (Fr,q(ap, a4) _ Frq(@k: bs) ) (14) 
q#k , 
where 
Pat x ,Y) 
by bea peta ba 

a a a 

ay, O41 


and the short form a#_, denotes the vector (21,..,U,—1,@41,--a)’ in (d — 1)- 
dimensions and x_,,, denotes the (d — 2)-dimensional vector (%1,...,%%-1, 

Lbet1, ++) Lq—1,Lq41; ++, La)’ for k £ g. The above equation (14) is deduced from 
Lee (1979, pp. 167). Note that for all t, = 0 the term Fy,,(x,y) will be the 
bivariate marginal density for which we will give a formula in the next section. 
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Subsequently, OF (Ph) can be obtained by substituting a; by b7,. From 12 — 


(15) at all t, =0,k ie 2,...,d, the second moment is 


_ F(t) 
~ “At;dt; 0 





E(X;X;) 


a 
; F; — 6, Fi,(b 
SoS ose (anFy (ax) — be FR (be)) 





k=1 kk 
g oO oO 
hack 
+I oD (og — EAH) [gles ag) ~ Fhgaes by) 
k=l  q#k kk 
F (Firqg( brs aq) - Frq(de: bq))] . (16) 


Having derived expressions for the first and second moments for double trun- 
cation in case of 4 = 0, we will now generalize to all uw. If Y ~ N(wu,™) with 
a* <y<b',thn X=Y—-—ypr~N(0,¥) witha=a*-p<ar<b-p=b 
and E(Y) = E(X) + pw and Cov(Y) = Cov(X). Equations 11 and 16 can then 
be used to compute E(X) and Cov(X). Hence, for general py, the first moment 
is 


d 
E(Y;) =) ois (Fi.(ax) — Fi(be)) + pie (17) 
k=l 


The covariance matrix 


is invariant to the shift in location. 

The equations 17 and 18 in combination with 11 and 16 form our derived 
result allow the calculation of the truncated mean and truncated variance for 
general double truncation. A formula for the term Fh 4(xx,%q), the bivariate 
marginal density, will be given in the next section. 

We have implemented the moment calculation for mean vector mean, co- 
variance matrix sigma and truncation vectors lower and upper as a function 
mtmvnorm(mean, sigma, lower, upper) in the R package tmvtnorm (Wilhelm 
& Manjunath, 2010, 2012), where the code is open source. In Section 5, we will 
show an example of this function. 


4 Bivariate Marginal Density Computation 


In order to compute the bivariate marginal density in this section, we follow 
Tallis (1961, p. 223) and Leppard and Tallis (1989) that implicitly used the 
bivariate marginal density as part of the moments calculation for single trunca- 
tion, evaluated at the integration bounds. However, we extend it to the doubly 
truncated case and state the function for all points within the support region. 
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Without loss of generality we use a z-transformation for all variates x = 
(a1,...,2q)/ as well as for all lower and upper truncation points a = (a 1,...,@a)! 
and b = (b1,..., ba)’, resulting in a N(0, R) distribution with correlation matrix 
R for the standardized untruncated variates. In this section we treat all variables 
as if they are z-transformed, leaving the notation unchanged. 

For computing the bivariate marginal density Fy,(%q,2,) with ag < xq < 
bg, Gr < tr < by,g #1, we use the fact that for truncated normal densities the 
conditional densities are also truncated normal. The following relationship holds 
for x, z5 € R¢~? conditionally on lq = Cg and z, =c, (s #qF 7): 





at pa(ts,Xq Cq, Ur cr; R) a1 p(€q; Cr} Par )Pa—2(2s3 Rar); (19) 





where 





25 = (5 — BsgrCq — Borger) / f= £24)(1 = Png) (20) 


and R,, is the matrix of the second-order partial correlation coefficients for 
SAGA. Bsqr and Bsrq are the partial regression coefficients of x, on #, and 
x, , respectively, and Ps,.q is the partial correlation coefficient between x, and 
x, for fixed xq. 


Integrating out (d — 2) variables x, leads to Fy,,(ag,%,) as a product of a 
bivariate normal density y(a,,z,) and a (d—2)-dimension normal integral ®g_2: 


bi bg-1 bg41 bp-1 
Fe (tg = Cle C) = / | oI 
ay @g—-1 JAg41 Qr—1 


br4i ba 
/ af PaR(2s,Cq, Cr) dX s 
Arti ad 


= a7 "—(Cq, Cr} Par )Pa-2( A; B4,; Ror) (21) 


where AZ, and B4, denote the lower and upper integration bounds of @g_2 given 
Lq = Cq and Lp = Cp: 





Al, = (0s — Boareq — Bor-atr)/ sf (1 — p2_)(1 - P2rq) (22) 


Ba, = (bs — Bogreg — Ber.atr)/ (1 — 024)(1 - Png): (23) 





The computation of Fy,-(#,2,) just needs the evaluation of the normal integral 
@a—2 in d — 2 dimensions, which is readily available in most statistics software 
packages, for example, as the function pmvnorm() in the R package mvtnorm 
(Genz et al., 2012). The bivariate marginal density function dtmvnorm(x, mean, 
sigma, lower, upper, margin=c(q,r)) is also part of the R package tmvtnorm 
(Wilhelm & Manjunath, 2010, 2012), where readers can find the source code as 
well as help files and additional examples. 
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5 Numerical Examples 


Example 1 


We will use the following bivariate example with = (0.5,0.5)’ and covariance 


matrix 
1 1.2 
E=(15 9) 


as well as lower and upper truncation points a = (—1,—oo)', b = (0.5, 1)’, ie. x1 
is doubly, while x2 is singly truncated. The bivariate marginal density F,,(z, y) 
is the density function itself and is shown in figure 1, and the one-dimensional 
densities Fi,(%) (k = 1,2) are shown in in figure 2. 


bivariate marginal density (x1,X2) 

















+ 
N A 
x of 

N 

] 

YH 
] ] ] T T 
-4 -2 0 2 4 

x4 


Figure 1. Contour plot for the bivariate truncated density function 


The moment calculation for our example can be performed in R as 


> library (tmvtnorm) 

> mu <- c(0.5, 0.5) 

> sigma <- matrix(c(1, 1.2, 1.2, 2), 2, 2) 
>a <- c(-1, -Inf) 

>b <- c(0.5, 1) 

> moments <- mtmvnorm(mean=mu, sigma=sigma, 
> lower=a, upper=b) 
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Figure 2. Marginal densities F(x) (k = 1,2) for v1 and x2 obtained from Kernel 
density estimation of random samples and from direct evaluation of F(x) 


which leads to the results w* = (—0.152, —0.388)’ and covariance matrix 


5 — (0-163 0.161 
~ (0.161 0.606 } ° 


The trace plots in figures 3 and 4 show the evolution of a Monte Carlo esti- 
mate for the elements of the mean vector and the covariance matrix respectively 
for growing sample sizes. Furthermore, the 95% confidence interval obtained 
from Monte Carlo using the full sample of 10000 items is shown. All confidence 
intervals contain the true theoretical value, but Monte Carlo estimates still show 
substantial variation even with a sample size of 10000. Simulation from a trun- 
cated multivariate distribution and calculating the sample mean or the sample 
covariance also leads to consistent estimates of u* and &*. Since the rate of 
convergence of the MC estimator is O(,/n), one has to ensure sufficient Monte 
Carlo iterations in order to have a good approximation or to choose variance 
reduction techniques. 


Example 2 
Let ps = (0,0,0)',the covariance matrix 


1.1 1.2 0 
Sr i a2 2 —0.8 
0 —0.8 3 
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Figure 3. Trace plots of the Monte Carlo estimator for p* 


and the lower and upper truncation points a = (—1, —co, —oo)’ and b = (0.5, 00, oo)’. 


Then the only truncated variable is x1, which is uncorrelated with x3. Our for- 
mula results in u* = c(—0.210, —0.229, 0)’ and 


0.174 0.190 0.0 
x* = | 0.190 0.898 —0.8 


0 -0.8 3.0 
For the special case of only k < d truncated variables (21,...,2,), the remaining 
d —k variables (a,41,...,a) can be regressed on the truncated variables, and 


a simple formula for the mean and covariance matrix can be given (see Johnson 


& Kotz, 1971, p. 70). 


Let the covariance matrix © of (#1,...,2q) be partitioned as 


Vit Vie 
> — 24 
ee ) - 
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Monte Carlo estimator for Cov(x;, x1) Monte Carlo estimator for Cov(x,, Xo) Monte Carlo estimator for Cov(xs, Xz) 
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Figure 4. Trace plots of the Monte Carlo estimator for the 3 elements of * (07, 


a” * * 
012 = 021 and 032) 


where V1; denotes the k x k covariance matrix of (#1,...,2,). The mean vector? 
and the covariance matrix &* of all d variables can be computed as 


(€,,€, Vi Viz) (25) 


and 





ye ( Un eee) (26) 
VaVi7Uu V22 — Val(Viy ViTUuVy))Vi2 


where €; and Uj; are the mean and covariance of the (#1,...,2%) after trunca- 
tion. 
The mean and standard deviation for the univariate truncated normal x are 


& =pt=oun Pur,011 (a1) _ Pur ors (Or) 
P yy or, (O1) — Pyy,ors (41) 
oft, =0114 oy eee (1) = b1 Pry ors (br) 
Py ,oy, (01) _ Py, ,01; (41) 








3 The formula for the truncated mean given in Johnson and Kotz (1971, p.70) is 
only valid for a zero-mean vector or after demeaning all variables appropriately. For 


non-zero means p= (fy, M2)’ it will be (€4, fy + (1 — #1) Vir Viz). 
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Letting Ui; = of, and inserting €; and Uj into equations 25 and 26, one 
can verify our formula and the results for w* and X*. However, the crux in 
using the Johnson/Kotz formula is the need to first compute the moments of the 
truncated variables (a1,...,2,) for k > 2. But this has been exactly the subject 
of our paper. 


6 Moment Calculation and Conditional Independence 


In this section we establish a link between our moment calculation and the the- 
ory of graphical models (Edwards, 1995; Lauritzen, 1996; Whittaker, 1990). We 
present some properties of the inverse covariance matrix and show how the de- 
pendence structure of variables is affected after selection. 


Graphical modeling uses graphical representations of variables as nodes in a 
graph and dependencies among them as edges. A key concept in graphical mod- 
eling is the conditional independence property. Two variables x and y are condi- 
tionally independent given a variable or a set of variables z (notation x IL y|z), 
when «x and y are independent after partialling out the effect of z. For condition- 
ally independent x and y, the edge between them in the graph is omitted and 
the joint density factorizes as f(x, y|z) = f(alz)f(ylz). 


Conditional independence is equivalent to having zero elements Qz, in the 
inverse covariance matrix Q = 7! as well as having a zero partial covari- 
ance/correlation between x and y given the remaining variables: 





a It y|Rest Qery 0 Pxy.Rest = 0. 





Both marginal independence and conditional independence between variables 
simplify the computations of the truncated covariance in equation 16. In the pres- 
ence of conditional independence of 7 and j given q, the terms 04; —ig0G, %qj = 0 
vanish as they reflect the partial covariance of i and 7 given q. 

As has been shown by Marchetti and Stanghellini (2008), the conditional 
independence property is preserved after selection, i.e. the inverse covariance 
matrices Q and * before and after truncation share the same zero-elements. 
We prove that many elements of the precision matrix are invariant to truncation. 
For the case of k < d truncated variables, we define the set of truncated variables 
with T = {x1,...,x,}, and the remaining d—k variables as S = {xp41,...,2a}- 
We can show that the off-diagonal elements 2; are invariant after truncation 
forte TUS andj €S: 


qq 


Proposition 1. The off-diagonal elements Q:,; and the diagonal elements Q;,; 
are invariant after truncation forte TUS andj eé S. 


Proof. The proof is a direct application of the Johnson/Kotz formula in equation 
26 in the previous section. As a result of the formula for partitioned inverse 
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matrices (Greene, 2003, section A.5.3) , the corresponding inverse covariance 
matrix 2 of the partitioned covariance matrix & is 


Q= Cee + Vi2F2V2iVi7) ee 


is 27 
—F2V21Vi; F, (27) 


with Fy = (V22 — VaVq Viz). 


Inverting the truncated covariance matrix }* in equation 26 using the for- 
mula for the partitioned inverse leads to the truncated precision matrix 


o* = Ce + Vi) Vi2F2VaiVit sae) (28) 


—F2V21Vi, F, 


where the 0}, and 93, elements are the same as Qy2 and Q2; respectively. The 
same is true for the elements in 35, especially the diagonal elements in Q5,. 














Here, we prove this invariance property only for a subset of truncated vari- 
ables. Based on our experiments we conjecture that the same is true also for the 
case of full truncation (i.e. all off-diagonal elements in Q},). However, we do not 
give a rigorous proof here and leave it to future research. 


Example 3 


We illustrate the invariance of the elements of the inverse covariance matrix with 
the famous mathematics marks example used in Whittaker (1990) and Edwards 
(1995, p.49). The independence graph of the five variables (W,V,X, Y,Z) in 
this example takes the form of a butterfly as shown in below. 


vectors (W) statistics (Z) 


algebra (X 


mechanics (V) analysis (Y) 


Here, we have the conditional independence (W,V) IL (Y,Z)|X. A corre- 
sponding precision matrix might look like (sample data; zero-elements marked 
as ”.”): 


1 02 03 
02: doa 
Q=|03-01 1040.5 (29) 
0.4 10.2 


0.50.2 1 
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After truncation in some variables (for example (W,V,X) as —2 < W <1, 
-1<V<1,0< X < 1), we apply equation 16 to compute the truncated 
second moment and the inverse covariance matrix as: 





1.88 0.2 0.3 
0.2 3.45 —0.1 
O* =] 0.3 -0.1 12.67 0.40.5 (30) 
0.4 10.2 
0.50.2 1 


The precision matrix * after selection differs from Q only in the diagonal 
elements of (W,V,X). From 2*, we can read how partial correlations between 
variables have changed due to the selection process. 

Each diagonal element 2,,, of the precision matrix is the inverse of the par- 
tial variance after regressing on all other variables (Whittaker, 1990, p. 143). 
Since only those diagonal elements in the precision matrix for the k < d of the 
truncated variables will change after selection, this leads to the idea to just com- 
pute these k elements after selection rather than the full k(& + 1)/2 symmetric 
elements in the truncated covariance matrix and applying the Johnson/Kotz 
formula for the remaining d— k variables. However, the inverse partial variance 
of a scalar y given the remaining variables X = {x,,...,xa}\y 





* ~ * *— 1 
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still requires the truncated covariance results derived in Section 3. 


7 Summary 


In this paper, we derived a formula for the first and second moments of the dou- 
bly truncated multivariate normal distribution and for their bivariate marginal 
density. An implementation for both formulas has been made available in the 
R statistics software as part of the tmvtnorm package. We linked our results to 
the theory of graphical models and proved an invariance property for elements 
of the precision matrix. Further research can deal with other types of truncation 
(e.g. elliptical). Another line of research can look at the moments of the dou- 
bly truncated multivariate Student-t distribution, which contains the truncated 
multivariate normal distribution as a special case. 
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